<!DOCTYPE html><html lang="zh-CN" data-theme="light"><head><meta charset="UTF-8"><meta http-equiv="X-UA-Compatible" content="IE=edge"><meta name="viewport" content="width=device-width,initial-scale=1"><title>弹簧质点系统与时间积分 | Tom's develop Blog</title><meta name="keywords" content="弹簧质点系统,时间积分"><meta name="author" content="Tom"><meta name="copyright" content="Tom"><meta name="format-detection" content="telephone=no"><meta name="theme-color" content="#ffffff"><meta http-equiv="Cache-Control" content="no-transform"><meta http-equiv="Cache-Control" content="no-siteapp"><meta name="description" content="弹簧质点系统与各种时间积分方法">
<meta property="og:type" content="article">
<meta property="og:title" content="弹簧质点系统与时间积分">
<meta property="og:url" content="https://tomsworkspace.github.io/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/index.html">
<meta property="og:site_name" content="Tom&#39;s develop Blog">
<meta property="og:description" content="弹簧质点系统与各种时间积分方法">
<meta property="og:locale" content="zh_CN">
<meta property="og:image" content="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/09/弹簧质点系统与时间积分/figure1.png">
<meta property="article:published_time" content="2021-08-09T13:39:12.000Z">
<meta property="article:modified_time" content="2021-09-13T12:20:27.206Z">
<meta property="article:author" content="Tom">
<meta property="article:tag" content="弹簧质点系统">
<meta property="article:tag" content="时间积分">
<meta name="twitter:card" content="summary">
<meta name="twitter:image" content="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/09/弹簧质点系统与时间积分/figure1.png"><link rel="shortcut icon" href="https://cdn.jsdelivr.net/gh/TOMsworkspace/BlogHelper@1.0/avatar.png"><link rel="canonical" href="https://tomsworkspace.github.io/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/"><link rel="preconnect" href="//cdn.jsdelivr.net"/><link rel="preconnect" href="//busuanzi.ibruce.info"/><link rel="preconnect" href="//zz.bdstatic.com"/><link rel="stylesheet" href="/css/index.css"><link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/@fortawesome/fontawesome-free/css/all.min.css"><link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/@fancyapps/fancybox@latest/dist/jquery.fancybox.min.css"><script>var GLOBAL_CONFIG = { 
  root: '/',
  algolia: undefined,
  localSearch: undefined,
  translate: {"defaultEncoding":2,"translateDelay":0,"msgToTraditionalChinese":"繁","msgToSimplifiedChinese":"簡"},
  noticeOutdate: undefined,
  highlight: {"plugin":"highlighjs","highlightCopy":true,"highlightLang":true},
  copy: {
    success: '复制成功',
    error: '复制错误',
    noSupport: '浏览器不支持'
  },
  relativeDate: {
    homepage: false,
    post: false
  },
  runtime: '天',
  date_suffix: {
    just: '刚刚',
    min: '分钟前',
    hour: '小时前',
    day: '天前',
    month: '个月前'
  },
  copyright: {"limitCount":50,"languages":{"author":"作者: Tom","link":"链接: ","source":"来源: Tom's develop Blog","info":"著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。"}},
  ClickShowText: {"text":"富强,民主,文明,和谐,自由,平等,公正,法治,爱国,敬业,诚信,友善","fontSize":"15px"},
  lightbox: 'mediumZoom',
  Snackbar: undefined,
  justifiedGallery: {
    js: 'https://cdn.jsdelivr.net/npm/justifiedGallery/dist/js/jquery.justifiedGallery.min.js',
    css: 'https://cdn.jsdelivr.net/npm/justifiedGallery/dist/css/justifiedGallery.min.css'
  },
  isPhotoFigcaption: false,
  islazyload: true,
  isanchor: false
};

var saveToLocal = {
  set: function setWithExpiry(key, value, ttl) {
    const now = new Date()
    const expiryDay = ttl * 86400000
    const item = {
      value: value,
      expiry: now.getTime() + expiryDay,
    }
    localStorage.setItem(key, JSON.stringify(item))
  },

  get: function getWithExpiry(key) {
    const itemStr = localStorage.getItem(key)

    if (!itemStr) {
      return undefined
    }
    const item = JSON.parse(itemStr)
    const now = new Date()

    if (now.getTime() > item.expiry) {
      localStorage.removeItem(key)
      return undefined
    }
    return item.value
  }
}</script><script id="config_change">var GLOBAL_CONFIG_SITE = { 
  isPost: true,
  isHome: false,
  isHighlightShrink: false,
  isToc: true,
  postUpdate: '2021-09-13 20:20:27'
}</script><noscript><style type="text/css">
  #nav {
    opacity: 1
  }
  .justified-gallery img {
    opacity: 1
  }

  #recent-posts time,
  #post-meta time {
    display: inline !important
  }
</style></noscript><script>(function () {  window.activateDarkMode = function () {
    document.documentElement.setAttribute('data-theme', 'dark')
    if (document.querySelector('meta[name="theme-color"]') !== null) {
      document.querySelector('meta[name="theme-color"]').setAttribute('content', '#0d0d0d')
    }
  }
  window.activateLightMode = function () {
    document.documentElement.setAttribute('data-theme', 'light')
   if (document.querySelector('meta[name="theme-color"]') !== null) {
      document.querySelector('meta[name="theme-color"]').setAttribute('content', '#ffffff')
    }
  }
  const autoChangeMode = 'true'
  const t = saveToLocal.get('theme')
  if (autoChangeMode === '1') {
    const isDarkMode = window.matchMedia('(prefers-color-scheme: dark)').matches
    const isLightMode = window.matchMedia('(prefers-color-scheme: light)').matches
    const isNotSpecified = window.matchMedia('(prefers-color-scheme: no-preference)').matches
    const hasNoSupport = !isDarkMode && !isLightMode && !isNotSpecified
    if (t === undefined) {
      if (isLightMode) activateLightMode()
      else if (isDarkMode) activateDarkMode()
      else if (isNotSpecified || hasNoSupport) {
        const now = new Date()
        const hour = now.getHours()
        const isNight = hour <= 6 || hour >= 18
        isNight ? activateDarkMode() : activateLightMode()
      }
      window.matchMedia('(prefers-color-scheme: dark)').addListener(function (e) {
        if (saveToLocal.get('theme') === undefined) {
          e.matches ? activateDarkMode() : activateLightMode()
        }
      })
    } else if (t === 'light') activateLightMode()
    else activateDarkMode()
  } else if (autoChangeMode === '2') {
    const now = new Date()
    const hour = now.getHours()
    const isNight = hour <= 6 || hour >= 18
    if (t === undefined) isNight ? activateDarkMode() : activateLightMode()
    else if (t === 'light') activateLightMode()
    else activateDarkMode()
  } else {
    if (t === 'dark') activateDarkMode()
    else if (t === 'light') activateLightMode()
  }})()</script><meta name="generator" content="Hexo 5.2.0"></head><body><div id="sidebar"><div id="menu-mask"></div><div id="sidebar-menus"><div class="author-avatar"><img class="avatar-img" data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/BlogHelper@1.0/avator.jpg" onerror="onerror=null;src='/img/friend_404.gif'" alt="avatar"/></div><div class="site-data"><div class="data-item is-center"><div class="data-item-link"><a href="/archives/"><div class="headline">文章</div><div class="length-num">38</div></a></div></div><div class="data-item is-center"><div class="data-item-link"><a href="/tags/"><div class="headline">标签</div><div class="length-num">73</div></a></div></div><div class="data-item is-center"><div class="data-item-link"><a href="/categories/"><div class="headline">分类</div><div class="length-num">12</div></a></div></div></div><hr/><div class="menus_items"><div class="menus_item"><a class="site-page" href="/"><i class="fa-fw fas fa-home"></i><span> Home</span></a></div><div class="menus_item"><a class="site-page" href="/archives/"><i class="fa-fw fas fa-archive"></i><span> Archives</span></a></div><div class="menus_item"><a class="site-page" href="/tags/"><i class="fa-fw fas fa-tags"></i><span> Tags</span></a></div><div class="menus_item"><a class="site-page" href="/categories/"><i class="fa-fw fas fa-folder-open"></i><span> Categories</span></a></div><div class="menus_item"><a class="site-page" href="/link/"><i class="fa-fw fas fa-link"></i><span> Link</span></a></div><div class="menus_item"><a class="site-page" href="/about/"><i class="fa-fw fas fa-heart"></i><span> About</span></a></div></div></div></div><div id="body-wrap"><header class="post-bg" id="page-header" style="background-image: url(https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/09/弹簧质点系统与时间积分/figure1.png)"><nav id="nav"><span id="blog_name"><a id="site-name" href="/">Tom's develop Blog</a></span><span id="menus"><div class="menus_items"><div class="menus_item"><a class="site-page" href="/"><i class="fa-fw fas fa-home"></i><span> Home</span></a></div><div class="menus_item"><a class="site-page" href="/archives/"><i class="fa-fw fas fa-archive"></i><span> Archives</span></a></div><div class="menus_item"><a class="site-page" href="/tags/"><i class="fa-fw fas fa-tags"></i><span> Tags</span></a></div><div class="menus_item"><a class="site-page" href="/categories/"><i class="fa-fw fas fa-folder-open"></i><span> Categories</span></a></div><div class="menus_item"><a class="site-page" href="/link/"><i class="fa-fw fas fa-link"></i><span> Link</span></a></div><div class="menus_item"><a class="site-page" href="/about/"><i class="fa-fw fas fa-heart"></i><span> About</span></a></div></div><span class="close" id="toggle-menu"><a class="site-page"><i class="fas fa-bars fa-fw"></i></a></span></span></nav><div id="post-info"><h1 class="post-title">弹簧质点系统与时间积分</h1><div id="post-meta"><div class="meta-firstline"><span class="post-meta-date"><i class="far fa-calendar-alt fa-fw post-meta-icon"></i><span class="post-meta-label">发表于</span><time class="post-meta-date-created" datetime="2021-08-09T13:39:12.000Z" title="发表于 2021-08-09 21:39:12">2021-08-09</time><span class="post-meta-separator">|</span><i class="fas fa-history fa-fw post-meta-icon"></i><span class="post-meta-label">更新于</span><time class="post-meta-date-updated" datetime="2021-09-13T12:20:27.206Z" title="更新于 2021-09-13 20:20:27">2021-09-13</time></span><span class="post-meta-categories"><span class="post-meta-separator">|</span><i class="fas fa-inbox fa-fw post-meta-icon"></i><a class="post-meta-categories" href="/categories/OpenSourceSummer2021/">OpenSourceSummer2021</a><i class="fas fa-angle-right post-meta-separator"></i><i class="fas fa-inbox fa-fw post-meta-icon"></i><a class="post-meta-categories" href="/categories/Computer-Graphics/">Computer Graphics</a><i class="fas fa-angle-right post-meta-separator"></i><i class="fas fa-inbox fa-fw post-meta-icon"></i><a class="post-meta-categories" href="/categories/Physics-engine/">Physics engine</a></span></div><div class="meta-secondline"> <span class="post-meta-separator">|</span><span class="post-meta-wordcount"><i class="far fa-file-word fa-fw post-meta-icon"></i><span class="post-meta-label">字数总计:</span><span class="word-count">1.6k</span><span class="post-meta-separator">|</span><i class="far fa-clock fa-fw post-meta-icon"></i><span class="post-meta-label">阅读时长:</span><span>7分钟</span></span><span class="post-meta-separator">|</span><span class="post-meta-pv-cv"><i class="far fa-eye fa-fw post-meta-icon"></i><span class="post-meta-label">阅读量:</span><span id="busuanzi_value_page_pv"></span></span><span class="post-meta-separator">|</span><span class="post-meta-commentcount"><i class="far fa-comments fa-fw post-meta-icon"></i><span class="post-meta-label">评论数:</span><a href="/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/#post-comment"><span class="gitalk-comment-count comment-count"></span></a></span></div></div></div></header><main class="layout" id="content-inner"><div id="post"><article class="post-content" id="article-container"><h2 id="弹簧质点系统"><a href="#弹簧质点系统" class="headerlink" title="弹簧质点系统"></a>弹簧质点系统</h2><p>一个模拟变形物体最简单的方法就是将其表示为弹簧质点系统（Mass Spring Systems）。一个弹簧质点包含了一系列由多个弹簧连接起来的质点，这样的系统的物理属性非常直接，模拟程序也很容易编写。</p>
<p>虽然模型简单，但是也带来了一些问题：</p>
<p>1.物体的行为依赖于弹簧系统的设置方法；</p>
<p>2.很难通过调整弹簧系数来得到想要的结果；</p>
<p>3.弹簧质点网格不能直接获取体效果。</p>
<p>在很多的应用中这些缺点可以忽略，在这种场合下，弹簧质点网格是最好的选择，因为够快够简单。弹簧质点系统可用于模拟绳索、布料、头发等弹性物体。</p>
<p>力是改变物体运动状态的原因，在这个系统中，主要有两种力，一是弹簧的弹力和阻尼力。</p>
<p>对于连接两个质点的一个弹簧，弹力是：<br><img src= "/img/loading.gif" data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/figure1.png" alt="弹力"></p>
<p>阻尼力可以这样算：<br><img src= "/img/loading.gif" data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/figure2.png" alt="阻尼力"></p>
<p>假设系统中有N个质点，质量为$m_i$，位置为$x_i$，速度为$v_i$ , $1 &lt; i &lt; N$.</p>
<p>这些质点由一组弹簧S连接，弹簧参数为（$i$, $j$, $l_0$, $k_s$, $k_d$）。$i$,$j$为连接的弹簧质点，$l_0$为弹簧完全放松时的长度，$k_s$为弹簧弹性系数，$k_d$为阻尼系数，由胡科定律知弹簧施加在两顶点上的力可以表示为：</p>
<p>$$ \overrightarrow{f}_i = \overrightarrow{f}^{s}(x_i,x_j)=k_s \frac{x_j - x_i}{|x_j - x _i|}(|x_j-x_i|-l_0) $$</p>
<p>$$ \overrightarrow{f}_j = \overrightarrow{f}^{s}(x_j,x_i)= -\overrightarrow{f}^{s}(x_i,x_j) = - \overrightarrow{f_i} $$</p>
<p>由受力守恒知$f_i+f_j = 0$. $f_i$和$f_j$的大小和弹簧伸长成正比关系。</p>
<p>对于阻尼的计算，除了与位移有关，还与质点速度有关：<br>$$ \overrightarrow{f_i} = \overrightarrow{f}^{d}(x_i,v_i,x_j,v_j)=k_d(v_j - v_i)\frac{x_j - x_i}{|x_j - x _i|} $$</p>
<p>$$ \overrightarrow{f_j} = \overrightarrow{f}^{s}(x_j,v_j,x_i,v_i) = -\overrightarrow{f_i} $$</p>
<p>大小和速度成正比，并且符合力守恒，则对于一个质点，其受合力方程为：</p>
<p>$$ \overrightarrow{f}(x_i,v_i)=\sum \overrightarrow{f}^{s}(x_i,x_j) + \sum \overrightarrow{f}^{d}(x_i,v_i,x_j,v_j) $$</p>
<p>这里 $j$ 为所有与质点 $i$ 存在弹簧连接的质点。后面讨论的运算都是矢量运算，为了方便就省略不写了。</p>
<p>在计算机模拟中，牛顿第二定律 $f = ma$ 是关键。在已知质量和外力的情况下，通过 $a=f/m $可以得到加速度，将二阶常微分方程写成两个一阶方程：<br>$$\frac{\mathrm{d}v}{dt} = \frac{f(x,v)}{m} $$</p>
<p>$$\frac{\mathrm{d}x}{dt} = v$$<br>可以得到解析解：</p>
<p>$$v(t) = v_0 + \int_{t_0}^{t}\frac{f(x,v)}{m} \mathrm{d}t$$</p>
<p>$$x(t) = x_0 + \int_{t_0}^{t}v(t) \mathrm{d}t$$</p>
<p>初始状态为 $v(t_0)$ = $v_0$, $x(t_0)=x_0$。积分将时间t内所有的变化加和，模拟的过程就是从$t_0$开始不断地计算$x(t)$和$v(t)$，然后更新质点的位置。</p>
<p>整个过程的伪代码如下：</p>
<figure class="highlight plain"><table><tr><td class="code"><pre><span class="line">&#x2F;&#x2F; initialization</span><br><span class="line"> forall particles i</span><br><span class="line">         initialize xi , vi and mi</span><br><span class="line"> endfor</span><br><span class="line">&#x2F;&#x2F; simulation loop</span><br><span class="line"> loop</span><br><span class="line">        forall particles i</span><br><span class="line">               fi ← fg + fcoll + ∑ f(xi , vi , x j , v j )</span><br><span class="line">        endfor</span><br><span class="line">        forall particles i</span><br><span class="line">             vi ← vi + ∆t fi &#x2F;mi</span><br><span class="line">             xi ← xi + ∆t vi</span><br><span class="line">         endfor</span><br><span class="line">         display the  every nth time</span><br><span class="line"> endloop</span><br></pre></td></tr></table></figure>

<h2 id="时间积分"><a href="#时间积分" class="headerlink" title="时间积分"></a>时间积分</h2><p>时间积分算法将这个积分的过程离散化，使用有限的步长去迭代下一时刻的状态。</p>
<h3 id="欧拉法"><a href="#欧拉法" class="headerlink" title="欧拉法"></a>欧拉法</h3><h4 id="前向欧拉方法-显式时间积分"><a href="#前向欧拉方法-显式时间积分" class="headerlink" title="前向欧拉方法 (显式时间积分)"></a>前向欧拉方法 (显式时间积分)</h4><p>$$v_{t + 1} = v_{t} + \Delta{t}\frac{f(x_t,v_t)}{m}$$</p>
<p>$$x_{t + 1} = x_{t} + \Delta{t}v_{t}$$</p>
<p>这个就是显式的欧拉解法，下一时刻的状态完全由当前状态决定。</p>
<h4 id="半隐式欧拉方法（显式时间积分）"><a href="#半隐式欧拉方法（显式时间积分）" class="headerlink" title="半隐式欧拉方法（显式时间积分）"></a>半隐式欧拉方法（显式时间积分）</h4><p>$$v_{t + 1} = v_{t} + \Delta{t}\frac{f(x_t,v_t)}{m}$$</p>
<p>$$x_{t + 1} = x_{t} + \Delta{t}v_{t + 1}$$</p>
<p>它和前向欧拉的差别很小。</p>
<h4 id="后向欧拉方法（隐式时间积分）"><a href="#后向欧拉方法（隐式时间积分）" class="headerlink" title="后向欧拉方法（隐式时间积分）"></a>后向欧拉方法（隐式时间积分）</h4><p>$$x_{t + 1} = x_{t} + \Delta{t}v_{t + 1}$$</p>
<p>$$v_{t + 1} = v_{t} + \Delta{t}M^{-1}f(x_{t + 1},v_{t + 1})$$</p>
<p>这里 $M$ 为一个 $3n * 3n$ 的对角矩阵，矩阵对角线上依次为各个质点的质量, $diag(M)=(m_1,m_1,m_1,m_2,m_2,m_2,…,m_n,m_n,m_n)$。</p>
<p>这是一个非线性的方程，先对 $f(x_{t + 1},v_{t + 1})$ 进行一阶泰勒展开，得</p>
<p>$$f(x_{t + 1},v_{t + 1}) = f(x_{t} + \Delta{x},v_{t} + \Delta{v})=f(x_{t},v_{t}) + \frac{\partial f}{\partial x}(x_t,v_t) \Delta{x} + \frac{\partial f}{\partial v}(x_t,v_t)\Delta{v} $$</p>
<p>由$\Delta{x} = x_{t+1}-x_{t} = \Delta{t}v_{t+1}$, $\Delta{v} = v_{t+1}-v_{t}$。得<br>$$f(x_{t + 1},v_{t + 1}) = f(x_{t},v_{t}) + \frac{\partial f}{\partial x}(x_t,v_t) \Delta{t} v_{t + 1} + \frac{\partial f}{\partial v}(x_t,v_t)(v_{t+1}-v_{t})$$</p>
<p>用 $I$ 表示单位矩阵，代入上面的式子，对上式进行移项，有</p>
<p>$$(I - \Delta{t}M^{-1}\frac{\partial f}{\partial v}(x_t,v_t) - (\Delta{t})^{2}M^{-1}\frac{\partial f}{\partial v}(x_{t},v_{t}))v_{t+1} = v_{t}+\Delta{t}M^{-1}(f(x_{t},v_{t})-\frac{\partial f}{\partial v}(x_{t},v_{t})v_{t})$$</p>
<p>$\frac{\partial f}{\partial v}$ 与 $\frac{\partial f}{\partial x}$ 的求法见<a target="_blank" rel="noopener" href="https://blog.mmacklin.com/2012/05/04/implicitsprings/">这里</a>。</p>
<p>上面的式子是一个形如，$AX = B$的方程组，$A$维度比较大，而且不一定可逆，一般使用<a target="_blank" rel="noopener" href="https://blog.csdn.net/weixin_40327927/article/details/88549879">雅克比迭代</a>等方法来解。</p>
<h3 id="中点法"><a href="#中点法" class="headerlink" title="中点法"></a>中点法</h3><p><a target="_blank" rel="noopener" href="https://en.wikipedia.org/wiki/Midpoint_method">中点法</a>是欧拉方法的改进，迭代方法如下：</p>
<p>$$x_{t + 1} = x_{t} + \Delta{t}v_{t + 1}$$</p>
<p>$$v_{t + 1} = v_{t} + \Delta{t}M^{-1}f(\frac{x_{t} + x_{t+1}}{2},\frac{v_{t} + v_{t+1}}{2})$$<br>求解与隐式欧拉相似。</p>
<h3 id="Heun法"><a href="#Heun法" class="headerlink" title="Heun法"></a>Heun法</h3><p><a target="_blank" rel="noopener" href="https://zh.wikipedia.org/wiki/Heun%E6%96%B9%E6%B3%95">Heun方法</a>是指改进或修改的显式欧拉方法，或类似的两阶段Runge-Kutta方法。</p>
<p>$$v_{t + 1} = v_{t} + \Delta{t}\frac{f(x_t,v_t)}{m}$$</p>
<p>$$x_{t + 1} = x_{t} + \frac{\Delta{t}}{2}(v_{t + 1} + v_{t})$$</p>
<h3 id="RK4法"><a href="#RK4法" class="headerlink" title="RK4法"></a>RK4法</h3><p>经典的<a target="_blank" rel="noopener" href="https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method">Runge Kutta方法</a>:</p>
<p>$$x_{t + 1} = x_{t} + \frac{1}{6}(k_{1}(x)+2k_2(x)+2k_3(x)+k_4(x))$$</p>
<p>$$v_{t + 1} = v_{t} + \frac{1}{6}(k_{1}(v)+2k_2(v)+2k_3(v)+k_4(v))$$</p>
<p>这里  </p>
<p>$$k_1(x) = \Delta{t}*v_t, k_1(v) = \Delta{t}*\frac{f(x_t,v_t)}{m}$$</p>
<p>$$k_2(x)=\Delta{t}*(v_t + \frac{k_1(v)}{2}), k_2(v) = \Delta{t}*\frac{f(x_t +\frac{k_1(x)}{2},v_t + \frac{k_1(v)}{2})}{m}$$</p>
<p>$$k_3(x)=\Delta{t}*(v_t + \frac{k_2(v)}{2}), k_3(v) = \Delta{t}*\frac{f(x_t +\frac{k_2(x)}{2},v_t + \frac{k_2(v)}{2})}{m}$$</p>
<p>$$k_4(x)=\Delta{t}*(v_t + k_3(v)), k_4(v) = \Delta{t}*\frac{f(x_t +k_3(x),v_t + k_3(v))}{m}$$</p>
<p><a target="_blank" rel="noopener" href="https://scicomp.stackexchange.com/questions/23929/equation-of-motion-by-rk4-method">代码实现</a></p>
<h3 id="Verlet法"><a href="#Verlet法" class="headerlink" title="Verlet法"></a>Verlet法</h3><p>基本的<a target="_blank" rel="noopener" href="https://en.wikipedia.org/wiki/Verlet_integration">Verlet积分</a></p>
<p>$$x_{t + 1} = 2x_{t} - x_{t-1}+ \Delta{t}v_{t + 1} + (\Delta t)^{2}\frac{f(x_t,v_t)}{m}$$</p>
<p>$$v_{t + 1} = \frac{x_{t+1}-x_{t}}{\Delta t} $$</p>
<p>参考<br>隐式欧拉：<br><a target="_blank" rel="noopener" href="https://blog.csdn.net/silangquan/article/details/12785001">https://blog.csdn.net/silangquan/article/details/12785001</a><br><a target="_blank" rel="noopener" href="https://zhuanlan.zhihu.com/p/148908332">https://zhuanlan.zhihu.com/p/148908332</a><br>RK4:<br><a target="_blank" rel="noopener" href="https://scicomp.stackexchange.com/questions/23929/equation-of-motion-by-rk4-method">https://scicomp.stackexchange.com/questions/23929/equation-of-motion-by-rk4-method</a><br>Jocabi迭代：<br><a target="_blank" rel="noopener" href="https://blog.csdn.net/weixin_40327927/article/details/88549879">https://blog.csdn.net/weixin_40327927/article/details/88549879</a></p>
</article><div class="post-copyright"><div class="post-copyright__author"><span class="post-copyright-meta">文章作者: </span><span class="post-copyright-info"><a href="mailto:undefined">Tom</a></span></div><div class="post-copyright__type"><span class="post-copyright-meta">文章链接: </span><span class="post-copyright-info"><a href="https://tomsworkspace.github.io/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/">https://tomsworkspace.github.io/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/</a></span></div><div class="post-copyright__notice"><span class="post-copyright-meta">版权声明: </span><span class="post-copyright-info">本博客所有文章除特别声明外，均采用 <a href="https://creativecommons.org/licenses/by-nc-sa/4.0/" target="_blank">CC BY-NC-SA 4.0</a> 许可协议。转载请注明来自 <a href="https://TOMsworkspace.github.io" target="_blank">Tom's develop Blog</a>！</span></div></div><div class="tag_share"><div class="post-meta__tag-list"><a class="post-meta__tags" href="/tags/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F/">弹簧质点系统</a><a class="post-meta__tags" href="/tags/%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/">时间积分</a></div><div class="post_share"><div class="social-share" data-image="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/09/弹簧质点系统与时间积分/figure1.png" data-sites="facebook,twitter,wechat,weibo,qq"></div><link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/social-share.js/dist/css/share.min.css"><script src="https://cdn.jsdelivr.net/npm/social-share.js/dist/js/social-share.min.js" defer></script></div></div><div class="post-reward"><div class="reward-button"><i class="fas fa-qrcode"></i> 打赏<div class="reward-main"><ul class="reward-all"><li class="reward-item"><a href="https://cdn.jsdelivr.net/gh/TOMsworkspace/BlogHelper@1.0/wp.jpg" target="_blank"><img class="post-qr-code-img" data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/BlogHelper@1.0/wp.jpg" alt="wechat"/></a><div class="post-qr-code-desc">wechat</div></li></ul></div></div></div><nav class="pagination-post" id="pagination"><div class="prev-post pull-left"><a href="/2021/09/11/%E5%BC%B9%E6%80%A7%E6%9C%89%E9%99%90%E5%85%83%E6%96%B9%E6%B3%95/"><img class="prev-cover" data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/09/11/弹性有限元方法/figure1.jpg" onerror="onerror=null;src='/img/404.jpg'"><div class="pagination-info"><div class="label">上一篇</div><div class="prev_info">弹性有限元方法</div></div></a></div><div class="next-post pull-right"><a href="/2021/08/08/OpenGL%E5%9D%90%E6%A0%87%E7%B3%BB%E7%BB%9F/"><img class="next-cover" data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/08/OpenGL坐标系统/figure1.png" onerror="onerror=null;src='/img/404.jpg'"><div class="pagination-info"><div class="label">下一篇</div><div class="next_info">OpenGL坐标系统</div></div></a></div></nav><hr/><div id="post-comment"><div class="comment-head"><div class="comment-headline"><i class="fas fa-comments fa-fw"></i><span> 评论</span></div></div><div class="comment-wrap"><div><div id="gitalk-container"></div></div></div></div></div><div class="aside_content" id="aside_content"><div class="card-widget card-info"><div class="card-content"><div class="card-info-avatar is-center"><img class="avatar-img" data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/BlogHelper@1.0/avator.jpg" onerror="this.onerror=null;this.src='/img/friend_404.gif'" alt="avatar"/><div class="author-info__name">Tom</div><div class="author-info__description">keep curiously!</div></div><div class="card-info-data"><div class="card-info-data-item is-center"><a href="/archives/"><div class="headline">文章</div><div class="length-num">38</div></a></div><div class="card-info-data-item is-center"><a href="/tags/"><div class="headline">标签</div><div class="length-num">73</div></a></div><div class="card-info-data-item is-center"><a href="/categories/"><div class="headline">分类</div><div class="length-num">12</div></a></div></div><a class="button--animated" id="card-info-btn" target="_blank" rel="noopener" href="https://github.com/TOMsworkspace"><i class="fab fa-github"></i><span>Follow Me</span></a><div class="card-info-social-icons is-center"><a class="social-icon" href="https://github.com/TOMsworkspace" target="_blank" title="Github"><i class="fab fa-github"></i></a><a class="social-icon" href="mailto:cqucscw@gmail.com" target="_blank" title="Email"><i class="fas fa-envelope"></i></a></div></div></div><div class="card-widget card-announcement"><div class="card-content"><div class="item-headline"><i class="fas fa-bullhorn card-announcement-animation"></i><span>公告</span></div><div class="announcement_content">This is my Blog</div></div></div><div class="sticky_layout"><div class="card-widget" id="card-toc"><div class="card-content"><div class="item-headline"><i class="fas fa-stream"></i><span>目录</span></div><div class="toc-content"><ol class="toc"><li class="toc-item toc-level-2"><a class="toc-link" href="#%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F"><span class="toc-number">1.</span> <span class="toc-text">弹簧质点系统</span></a></li><li class="toc-item toc-level-2"><a class="toc-link" href="#%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86"><span class="toc-number">2.</span> <span class="toc-text">时间积分</span></a><ol class="toc-child"><li class="toc-item toc-level-3"><a class="toc-link" href="#%E6%AC%A7%E6%8B%89%E6%B3%95"><span class="toc-number">2.1.</span> <span class="toc-text">欧拉法</span></a><ol class="toc-child"><li class="toc-item toc-level-4"><a class="toc-link" href="#%E5%89%8D%E5%90%91%E6%AC%A7%E6%8B%89%E6%96%B9%E6%B3%95-%E6%98%BE%E5%BC%8F%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86"><span class="toc-number">2.1.1.</span> <span class="toc-text">前向欧拉方法 (显式时间积分)</span></a></li><li class="toc-item toc-level-4"><a class="toc-link" href="#%E5%8D%8A%E9%9A%90%E5%BC%8F%E6%AC%A7%E6%8B%89%E6%96%B9%E6%B3%95%EF%BC%88%E6%98%BE%E5%BC%8F%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86%EF%BC%89"><span class="toc-number">2.1.2.</span> <span class="toc-text">半隐式欧拉方法（显式时间积分）</span></a></li><li class="toc-item toc-level-4"><a class="toc-link" href="#%E5%90%8E%E5%90%91%E6%AC%A7%E6%8B%89%E6%96%B9%E6%B3%95%EF%BC%88%E9%9A%90%E5%BC%8F%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86%EF%BC%89"><span class="toc-number">2.1.3.</span> <span class="toc-text">后向欧拉方法（隐式时间积分）</span></a></li></ol></li><li class="toc-item toc-level-3"><a class="toc-link" href="#%E4%B8%AD%E7%82%B9%E6%B3%95"><span class="toc-number">2.2.</span> <span class="toc-text">中点法</span></a></li><li class="toc-item toc-level-3"><a class="toc-link" href="#Heun%E6%B3%95"><span class="toc-number">2.3.</span> <span class="toc-text">Heun法</span></a></li><li class="toc-item toc-level-3"><a class="toc-link" href="#RK4%E6%B3%95"><span class="toc-number">2.4.</span> <span class="toc-text">RK4法</span></a></li><li class="toc-item toc-level-3"><a class="toc-link" href="#Verlet%E6%B3%95"><span class="toc-number">2.5.</span> <span class="toc-text">Verlet法</span></a></li></ol></li></ol></div></div></div><div class="card-widget card-recent-post"><div class="card-content"><div class="item-headline"><i class="fas fa-history"></i><span>最新文章</span></div><div class="aside-list"><div class="aside-list-item"><a class="thumbnail" href="/2021/09/27/%E6%B5%81%E4%BD%93%E6%A8%A1%E6%8B%9F%E4%B8%8ESPH%E6%96%B9%E6%B3%95/" title="流体模拟与SPH方法"><img data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/09/27/流体模拟与SPH方法/figure1.jpg" onerror="this.onerror=null;this.src='/img/404.jpg'" alt="流体模拟与SPH方法"/></a><div class="content"><a class="title" href="/2021/09/27/%E6%B5%81%E4%BD%93%E6%A8%A1%E6%8B%9F%E4%B8%8ESPH%E6%96%B9%E6%B3%95/" title="流体模拟与SPH方法">流体模拟与SPH方法</a><time datetime="2021-09-27T08:46:17.000Z" title="发表于 2021-09-27 16:46:17">2021-09-27</time></div></div><div class="aside-list-item"><a class="thumbnail" href="/2021/09/11/%E5%BC%B9%E6%80%A7%E6%9C%89%E9%99%90%E5%85%83%E6%96%B9%E6%B3%95/" title="弹性有限元方法"><img data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/09/11/弹性有限元方法/figure1.jpg" onerror="this.onerror=null;this.src='/img/404.jpg'" alt="弹性有限元方法"/></a><div class="content"><a class="title" href="/2021/09/11/%E5%BC%B9%E6%80%A7%E6%9C%89%E9%99%90%E5%85%83%E6%96%B9%E6%B3%95/" title="弹性有限元方法">弹性有限元方法</a><time datetime="2021-09-11T09:05:34.000Z" title="发表于 2021-09-11 17:05:34">2021-09-11</time></div></div><div class="aside-list-item"><a class="thumbnail" href="/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/" title="弹簧质点系统与时间积分"><img data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/09/弹簧质点系统与时间积分/figure1.png" onerror="this.onerror=null;this.src='/img/404.jpg'" alt="弹簧质点系统与时间积分"/></a><div class="content"><a class="title" href="/2021/08/09/%E5%BC%B9%E7%B0%A7%E8%B4%A8%E7%82%B9%E7%B3%BB%E7%BB%9F%E4%B8%8E%E6%97%B6%E9%97%B4%E7%A7%AF%E5%88%86/" title="弹簧质点系统与时间积分">弹簧质点系统与时间积分</a><time datetime="2021-08-09T13:39:12.000Z" title="发表于 2021-08-09 21:39:12">2021-08-09</time></div></div><div class="aside-list-item"><a class="thumbnail" href="/2021/08/08/OpenGL%E5%9D%90%E6%A0%87%E7%B3%BB%E7%BB%9F/" title="OpenGL坐标系统"><img data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/08/OpenGL坐标系统/figure1.png" onerror="this.onerror=null;this.src='/img/404.jpg'" alt="OpenGL坐标系统"/></a><div class="content"><a class="title" href="/2021/08/08/OpenGL%E5%9D%90%E6%A0%87%E7%B3%BB%E7%BB%9F/" title="OpenGL坐标系统">OpenGL坐标系统</a><time datetime="2021-08-08T12:57:28.000Z" title="发表于 2021-08-08 20:57:28">2021-08-08</time></div></div><div class="aside-list-item"><a class="thumbnail" href="/2021/08/06/Doxygen%E8%87%AA%E5%8A%A8%E7%94%9F%E6%88%90%E6%96%87%E6%A1%A3/" title="Doxygen自动生成文档"><img data-lazy-src="https://cdn.jsdelivr.net/gh/TOMsworkspace/TOMsworkspace.github.io/2021/08/06/Doxygen自动生成文档/Doxygen.png" onerror="this.onerror=null;this.src='/img/404.jpg'" alt="Doxygen自动生成文档"/></a><div class="content"><a class="title" href="/2021/08/06/Doxygen%E8%87%AA%E5%8A%A8%E7%94%9F%E6%88%90%E6%96%87%E6%A1%A3/" title="Doxygen自动生成文档">Doxygen自动生成文档</a><time datetime="2021-08-06T13:24:50.000Z" title="发表于 2021-08-06 21:24:50">2021-08-06</time></div></div></div></div></div></div></div></main><footer id="footer"><div id="footer-wrap"><div class="copyright">&copy;2019 - 2021 By Tom</div><div class="framework-info"><span>框架 </span><a target="_blank" rel="noopener" href="https://hexo.io">Hexo</a><span class="footer-separator">|</span><span>主题 </span><a target="_blank" rel="noopener" href="https://github.com/jerryc127/hexo-theme-butterfly">Butterfly</a></div><div class="footer_custom_text">Hi, welcome to my <a href="https://TOMsworkspace.github.io/">blog</a>!</div></div></footer></div><div id="rightside"><div id="rightside-config-hide"><button id="readmode" type="button" title="阅读模式"><i class="fas fa-book-open"></i></button><button id="translateLink" type="button" title="简繁转换">简</button><button id="darkmode" type="button" title="浅色和深色模式转换"><i class="fas fa-adjust"></i></button></div><div id="rightside-config-show"><button id="rightside_config" type="button" title="设置"><i class="fas fa-cog"></i></button><button id="chat_btn" type="button" title="rightside.chat_btn"><i class="fas fa-sms"></i></button><a id="to_comment" href="#post-comment" title="直达评论"><i class="fas fa-comments"></i></a><button id="go-up" type="button" title="回到顶部"><i class="fas fa-arrow-up"></i></button></div></div><div><script src="https://cdn.jsdelivr.net/npm/jquery@latest/dist/jquery.min.js"></script><script src="/js/utils.js"></script><script src="/js/main.js"></script><script src="/js/tw_cn.js"></script><script src="https://cdn.jsdelivr.net/npm/medium-zoom/dist/medium-zoom.min.js"></script><script src="https://cdn.jsdelivr.net/npm/instant.page/instantpage.min.js" type="module" defer></script><script src="https://cdn.jsdelivr.net/npm/vanilla-lazyload/dist/lazyload.iife.min.js"></script><div class="js-pjax"><script>if (!window.MathJax) {
  window.MathJax = {
    loader: {
      source: {
        '[tex]/amsCd': '[tex]/amscd'
      }
    },
    tex: {
      inlineMath: [ ['$','$'], ["\\(","\\)"]],
      tags: 'ams'
    },
    options: {
      renderActions: {
        findScript: [10, doc => {
          for (const node of document.querySelectorAll('script[type^="math/tex"]')) {
            const display = !!node.type.match(/; *mode=display/)
            const math = new doc.options.MathItem(node.textContent, doc.inputJax[0], display)
            const text = document.createTextNode('')
            node.parentNode.replaceChild(text, node)
            math.start = {node: text, delim: '', n: 0}
            math.end = {node: text, delim: '', n: 0}
            doc.math.push(math)
          }
        }, ''],
        addClass: [200,() => {
          document.querySelectorAll('mjx-container:not([display=\'true\']').forEach( node => {
            const target = node.parentNode
            if (!target.classList.contains('has-jax')) {
              target.classList.add('mathjax-overflow')
            }
          })
        }, '', false]
      }
    }
  }
  
  const script = document.createElement('script')
  script.src = 'https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js'
  script.id = 'MathJax-script'
  script.async = true
  document.head.appendChild(script)
} else {
  MathJax.startup.document.state(0)
  MathJax.texReset()
  MathJax.typeset()
}</script><link rel="stylesheet" type="text/css" href="https://cdn.jsdelivr.net/npm/katex@latest/dist/katex.min.css"><script src="https://cdn.jsdelivr.net/npm/katex-copytex@latest/dist/katex-copytex.min.js"></script><link rel="stylesheet" type="text/css" href="https://cdn.jsdelivr.net/npm/katex-copytex@latest/dist/katex-copytex.min.css"><script>$(function () {
  $('span.katex-display').wrap('<div class="katex-wrap"></div>')
})</script><script>function addGitalkSource () {
  const ele = document.createElement('link')
  ele.rel = 'stylesheet'
  ele.href= 'https://cdn.jsdelivr.net/npm/gitalk/dist/gitalk.min.css'
  document.getElementsByTagName('head')[0].appendChild(ele)
}

function loadGitalk () {
  function initGitalk () {
    var gitalk = new Gitalk({
      clientID: 'e50742231659674c31e3',
      clientSecret: '59038139e289bc384e484d341205bbb95a123209',
      repo: 'TOMsworkspace.github.io',
      owner: 'TOMsworkspace',
      admin: ['TOMsworkspace'],
      id: '3c18b8570748188c0020063fc88f6c0f',
      language: '-zh-CN# en, zh-CN, zh-TW, es-ES, fr, ru',
      perPage: 10,
      distractionFreeMode: false,
      pagerDirection: 'last',
      createIssueManually: true,
      updateCountCallback: commentCount
    })
    gitalk.render('gitalk-container')
  }

  if (typeof Gitalk === 'function') initGitalk()
  else {
    addGitalkSource()
    $.getScript('https://cdn.jsdelivr.net/npm/gitalk@latest/dist/gitalk.min.js', initGitalk)
  }
}

function commentCount(n){
  let isCommentCount = document.querySelector('#post-meta .gitalk-comment-count')
  if (isCommentCount) {
    isCommentCount.innerHTML= n
  }
}

if ('Gitalk' === 'Gitalk' || !false) {
  if (false) btf.loadComment(document.getElementById('gitalk-container'), loadGitalk)
  else loadGitalk()
} else {
  function loadOtherComment () {
    loadGitalk()
  }
}</script><script async src="//busuanzi.ibruce.info/busuanzi/2.3/busuanzi.pure.mini.js"></script></div><script defer="defer" id="fluttering_ribbon" mobile="false" src="https://cdn.jsdelivr.net/npm/butterfly-extsrc@1/dist/canvas-fluttering-ribbon.min.js"></script><script id="click-show-text" src="https://cdn.jsdelivr.net/npm/butterfly-extsrc@1/dist/click-show-text.min.js" async="async" mobile="false"></script><script>(function(){
  const bp = document.createElement('script');
  const curProtocol = window.location.protocol.split(':')[0];
  if (curProtocol === 'https'){
  bp.src = 'https://zz.bdstatic.com/linksubmit/push.js';
  }
  else{
  bp.src = 'http://push.zhanzhang.baidu.com/push.js';
  }
  bp.dataset.pjax = ''
  const s = document.getElementsByTagName("script")[0];
  s.parentNode.insertBefore(bp, s);
})()</script></div></body></html>